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A  KALMAN  FILTER  FOR  A  POISSON  SERIES  WITH  COVARIATES 
AND  LAPLACE-APPROXIMATION  INTEGRATION 

D.  P.  Gaver 
P.  A.  Jacobs 

0.     INTRODUCTION 

The  Poisson  model  is  an  initial  idealized,  but  plausible,  off-the-shelf  tool 
for  representing  point-process  data  of  nearly  any  kind,  cf.  Feller  (1966)  and 
Cox  and  Lewis  (1968).  However,  to  be  more  descriptive,  and  even  predictive, 
representation  of  non-homogeneity  in  space  or  time  may  be  needed.  For 
instance  the  occurrence  of  rare  events  in  space,  such  as  the  occurrence  of 
extreme  heights,  i.e.  above  a  level  of  Arctic  ice  along  a  transsect  or  encounters 
with  ice  keels  along  a  submarine  track  at  constant  (deep)  depths,  may  well 
appear  roughly  Poisson.  A  better  description  of  these  events  may  require 
more  detail  than  a  simple  mean  or  rate:  some  account  of  regional  and 
seasonal  variation  could  be  needed  for  true  accuracy;  for  instance  the 
intervention  of  natural  gaps  along  a  path  in  Arctic  ice  will  occur  if  the  latter 
crosses  leads  (open  water  amidst  an  ice  pack).  For  another  example,  demands 
for  spare  parts  in  a  logistics  system  often  are  roughly  Poisson,  or  compound 
Poisson,  but  with  a  mean  or  rate  that  changes  erratically  but  slowly  in  time. 
Demands  for  communication  and  computer  facilities  exhibit  a  similar 
temporal  pattern,  and  there  are  a  great  many  other  examples.  To  summarize, 
variation  in  a  fundamental  Poisson  rate  or  mean  is  very  often  encountered  in 
practice;  it  is  possible  that  if  this  variation  is  slow-moving  or  persistent 
enough  it  can  be  exploited  for  short-term  forecasting. 


In  this  paper  we  study  a  model-based  procedure  for  forecasting  in  the 
kind  of  environments  described.  The  model  introduced  allows  the  mean  or 
rate  of  the  Poisson  process  to  be  itself  a  random  process;  the  exponential  of  an 
AR/1  autoregressive  process.  In  addition  the  rate  is  influenced  by  a 
covariate.  We  then  recursively  update  the  parameter  estimates  using  an 
approximation  based  on  the  Laplace  method;  cf.  de  Bruijn  (1958).  The  approach 
resembles  that  of  Delampady,  Yee  and  Zidek  (1991),  but  frankly  heuristic 
methods  are  used  to  estimate  certain  of  the  underlying  parameters.  The 
methodology  is  checked  against  simulated  data  with  encouraging  results. 

1.     FORMULATION 

Consider  the  following  Poissonian  model  for  count  data: 

P{Y,  =ytkt,xt,frht,y0,...,yt_1}  =  exp{-htex*+>ity- — l—        (1.1) 

Vt  • 

where  •  > 

\h  =a\Lt-i+®t  (1-2) 

with  {©,}  independent  normal/Gaussian  random  variables  with  mean  0  and 
variance  W t  independent  of  {u.,-;  i  <  f-1},  {Y,;  i  <  t-1),  {xt},  {ht}  and  p.   Another 
hierarchical  time  series  model  for  count  data  can  be  found  in  Harvey  et  al. 
(1989). 

The  purpose  of  this  paper  is  to  suggest  a  Kalman  filter-like  procedure  to 
produce  successive  estimates  of  P  and  \it  as  new  data  becomes  available.  The 
procedure  is  based  on  a  Laplace  approximation  to  an  integral.  A  Bayesian 
approach  to  a  similar  problem  is  being  investigated  by  Delampaday  et  al. 


(1991).   A  Bayesian  approach  to  time  series  can  be  found  in  West  et  al.  (1989) 
and  for  a  more  recent  computational  approach  in  Carlin  et  al.  (1991). 

2.     AN  APPROXIMATE  UPDATING  PROCEDURE 

Assume  that  the  posterior  distribution  of  (p,  jifl)  given  {y„  i<  t-1)  is 
bivariate  normal  with  mean  (bt_v  m^),  Var[fV]  =  tt_v  Var[p.fl]  =  Ct-1  and 
Corr[p,uM]  =  pt_r 

Since  it  is  known  that 

the  prior  distribution  of  (P,jif)  is  bivariate  normal  with  mean  {bt_y  am,^), 
Var[p]  =  xt_v  Var[jit]  =  Rt  =  a  CM  +  Wt  and  Corr[p,jif]  =  rt  =  apf_1  y]Ct-i/Rt. 

The  forecast/prediction  distribution  of  Yf  in  terms  of  data  up  to  t-1  and 
the  covariate  value  at  t  is 

P{\t  =  yf|Y,  =  yiti<t-l,xifi  <  t,hi,i<t} 


=  jexp{.e{Xib+%} 


\xtb+z\. 


Vl 


*exp 


4M) 


-l 


{b-btf     2rt{b-b?)(z-m{)     (z-m{\ 
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VvTV^" 
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where  bf  =  bt.\  and mf  =  amt _j. 


We  now  approximate  the  integral  by  the  Laplace  method;  cf.  Easton 
(1991),  Cox  and  Hinkley  (1974),  de  Bruijn  (1958).  Let  the  exponent  of  the 
integrand  be 


g{b,z)  =  -  e^b+%  +  yt[xtb  *  z\  +  yf  In*,  -  A(l  -  rf2) " ! 
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where  X  is  a  constant. 
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Use  a  Newton  procedure  to  solve  the  system  of  equations 
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dz 
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for  mf  and  bt-  Solve  for  t,,  Q  and  pt  using  the  relation 
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The  posterior  distribution  for  (P,|if)  given  yt,  D t-\  is  approximated  by  a 
bivariate  normal  distribution  with  mean  (bt,  nit)  and  variance-covariance 
matrix 

Tf  Pt^cth 

PtJCt*t  ct 


*t- 


(2.10) 


2.1    Summary  of  the  Newton  Procedure 

0.  Start  with  estimates  of  the  parameters  of  the  prior  bivariate  distribution 
of  (P,jif);  that  is,  estimates  of  the  mean  (bt-\,  amt-i),  Corr(P,u.t)  =  rt,  Var(p)  =  rt-i 
Var[u.(]  =  Rt.  Set  bt  -  bM,  m  t  =  cuniv 

1.  Solve  the  system  of  linear  equations  (in  b  and  m) 

2  2 

d     /,o      n*  .    d        /,o      n\r,     ,.m  .   d 


„0 
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dzdb 


for  b,  m . 


2.  If  max  (fy  -b  \/bj,umt  -  m  /  m)  <  0.001,  set  mf  =  m  and  bt=  b  and  go  to 
3.  Otherwise  set  bt  =b  and  mt  =m  and  return  to  Step  1  unless  Step  1  has  been 
returned  to  49  times  in  which  case  set  mt=  mt  and  bt  =  bt;  go  to  3. 

3.  Return  (bt,tnt)  as  the  estimate  of  the  posterior  mean  of  the  bivariate 
normal  distribution  of  (P,[i/). 


4.  To  find  the  variance-covariance  matrix  of  the  posterior  distribution  of 
(p,uf)  evaluate 

^'^        -Skg{bt'mt) 

dzdb  dz 

set  it  equal  to 


and  solve  for  xt ,  pt  and  C/. 

2.2    Summary  of  Kalman  Procedure 

In  summary  the  approximate  Kalman  procedure  is  as  follows 

0.  Start  with  the  parameters  of  the  (approximate)  posterior  bivariate 
normal  distribution  of  ($,\Lt-\)  having  mean  (bt-\,  mt-i),  Corr((J,[if_i)  =  pt-\, 
Var[p]  =  tt-i  Var[jiM]  =  CM- 

1.  Update  Corr(P,p.f_i)  and  Var[p.f.i]  using  (1.2)  to  obtain  the  prior 
bivariate  normal  distribution  of  ($,\Lt)  having  mean  (bt-\,  amt-\),  Var[f$]  =  Xt-\, 
Var[jif]  ■  Rt  =  a  CM  +  Wtr  Corr(P,u.,)  ■  rt  =  apt-\  \lc^/Rt. 

2.  Observe  the  Poisson  count  yt. 

3.  Invoke  the  Newton  procedure  of  Section  2.1  to  obtain  estimates  of  the 
parameters  of  the  approximate  posterior  distribution  of  (P,|it)  given  past 
observations  and  the  new  observation  yt- 


k 


To  obtain  moments  for  the  Poisson  mean  Xt  -  hfixp{xfi+\it)  note  that  since 
the  distribution  of  (fit,  \it)  is  being  approximated  by  a  bivariate  normal 
distribution,  the  posterior  moments  of  Xt  are  approximately 

E[A*]  *  expjx^  +mf  +i[x?T,  +  Ct  +  2xtPtJ^Jq\ 

Var[X,]  *  /z2exp{2(x,&,  +  mf)  +  [x2t,  +Ct  +2xtptJ^Jc~t\} 
x[exp{A:2Tf  +  C,  +  2xtPt  yfcjc;}  - 1]. 

The  estimate  of  X,is  Xt  =  ebt%l  +m%. 

2.3   A  Simulation  Example 

In  the  example  P  =  0.5,  {xt}  takes  the  values  {0.25,  .5, 1}  over  and  over;  ht  -  1. 
{(Ot)  are  iid  normal  with  mean  0  and  variance  0.25  and  a  =  0.5.  Given  P  and  \Ltl 
Y(  has  a  Poisson  distribution  with  mean  Xt  =  ex^  ^f .   The  simulation  starts 

with  Ho  drawn  from  the  stationary  distribution  of  {\it}  a  normal  distribution 

2       1 
with  mean  0  and  variance  W/(l-a  )  =  ~.  Successive  u*  are  computed  as 

\it  =a\Lt.1+(at'/ 

the  Xt  is  computed  and  a  Poisson  random  number  with  mean  Xt  is  then 
generated.  The  simulation  data  are  the  Poisson  counts  and  {xt};  t  =  1,  ...,  100. 
The  random  numbers  were  generated  using  LLRANDOM  II;  cf.  Lewis  et  al. 
(1981). 

At  time  0,  the  Newton  procedure  is  initialized  at  fio  =  0,  j3  =  0,  po  =  tq  =  0, 
a  =  0.5,  W  =  0.25,  C0  =  0,  r\  =  0.25,  To  =  1;  note  the  a  and  W  are  assumed 
known.  The  data  point  yi  is  observed  and  the  Newton  procedure  is  used  to 
find  the  posterior  moments  [m\,b\,  p\,C\,  T\]. 


The  da  a  point  1/2  is  the  observed  and  the  Newton  procedure  is  started 
with  [m\,  b  1,  t\,  C\  +  7,  Ti]  and  used  to  find  [mi,  h,  Pi,  C2,  T2],  etc. 

Results  of  the  simulation  appear  in  Figures  1-3.  In  Figure  1  the  count  data 
yt  appears  along  with  the  true  \t  (dotted  line)  and  the  estimated  Xt  (solid  line). 
In  Figure  2  the  true  value  of  \it  (dotted  line)  and  estimated  value  of  \i{  (solid 
line)  appear.  The  estimated  values  of  the  standard  deviation  of  \Lt,  y[Ct ,  also 
appear  in  Figure  2.  The  estimated  value  of  /3,  P's  estimated  standard 
deviation  yjzt,  and  the  estimated  value  of  pt  appear  in  Figure  3.  Figures  1-3 
indicate  that  the  procedure  performs  satisfactorily.  The  apparent  oscillation 
in  some  of  the  figures,  (particularly  fit),  is  due  to  the  cyclic  nature  of  {xt}.  Not 
surprisingly,  the  estimated  values  of  X.t  and  \it  are  less  variable  than  the  true 
values.  They  appear  to  be  practically  acceptable. 

3.  APPROXIMATE  KALMAN  PROCEDURES  WHICH  INCORPORATE 
ESTIMATES  OF  THE  AUTOREGRESSIVE  PARAMETERS  a  AND  W: 
NAIVE  MOMENT  ESTIMATORS 

In  this  section  we  will  assume  that  {coj  are  independent  identically 

distributed  normal/Gaussian  random  variables  with  mean  0  and  variance  W 

with  co,  independent  of  {\L{;  i  <  M},  {Y-;  i  <  t -  1},  [ht]  and  p. 

3.1    Estimates  of  the  Autoregressive  Parameters  a  and  W 

If  {\it;  t  <  T)  were  observable,  then  one  could  estimate  a  and  W  by 
maximum  likelihood;  that  is,  the  likelihood  function  is 

L(W'°° = n^Jw  exp[  i(M(  •  aMt-i)2^}         (ai) 

or 


Now 


gives 


and 


from  which 


lnL(W,a)  -  l(W,a)  =  -£  (/x,  -  «M^-i)2^- 

T  T 

-— InW  --ln27T.  (3.2) 

2  2  *"      ; 


r\l  1 

—  =  -£fo  -aMf.i)Mf-i-  =  0  (3.3) 


T  T 

a(T)  =  5>t/^.i/5>?-i  (3.4) 


a/     i  i  ^  ,2    t  i 


fife  -«Mf-l)    -TiT7=0  (3-5) 


aw      2wzPi  2W 


7 

VV(T)  =  -£(/x,  -a(T)Mf.!)2.  (3.6) 

1  t=l 

Unfortunately  nfisnot  observable.  One  possible  estimate  of  \it  is  the 
posterior  mean  mt  of  subsections  (2.1)  and  (2.2).  In  this  case  the  corresponding 
estimators  of  a  and  W  are 

T  T 

«m(^)  =  I^^-l/I^2-l  (3-7) 

1    T  2 

Wm(T)=-^(mt-am(T)mt.1)Z.  (3.8) 

J  f=l 


Another  estimate  of  p.,  is 


My)  -  in 


1 
Vt  +- 


xtfy 


(3.9) 


where  bt  is  the  mean  of  the  approximate  posterior  distribution  of  {$  given 
{ys;  s  <  t}.  The  corresponding  estimates  of  a  and  W  are 

T  T 


and 


ayCn  -X(&(y)&-i(y))/X/**(y)' 


T 

^(T)4lfe(y)-«y(^-i(y))' 


(3.10) 


(3.11) 


The  estimates  (3.7)  -  (3.11)  can  be  recomputed  at  every  time  T.  Other 
similar  estimates  can  be  obtained  by  not  recomputing  the  estimates  at  every 
time  T;  for  example,  choose  an  integer  8  >  1,  put 

a(t;S)  =  a{{n-l)S;8)  (3.12) 


and 


if  (n  -1)5  <  f  <nS; 


W(t;8)  =W((n  -1)8,8) 


(3.13) 


n8 


n8 


d(n8;8)  =  £/2f&.i/£&-i 

1   "5   .       „  .2 

W(n8;8)  =  —  £  fof  -  a(n5;5)/if  .j) 


n«5 


(3.14) 


(3.15) 


f=l 


where  flt  is  an  estimate  of  \ir  Another  possibility  is  to  use  a  window  of  times 
to  compute  estimates  of  a  and  W. 
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The  following  is  a  summary  of  the  Kalman  procedure  of  Section  2  with 
the  addition  of  estimation  of  a  and  W. 

3.2  Summary  of  the  Kalman  Procedure  with  Estimation  of  a  and  W 

0.  Start  with  the  parameters  of  the  (approximate)  posterior  bivariate 
normal  distribution  of  (P,  \it_^)  having  mean  {bt_v  mt_i),  Corr(P,u.t_j)  =  pt_ 
v  VarfP]  =  xt_y  Var[nM]  =  Ct_r 

1.  Update  Corr(P,  |tt_1)  and  Var[p.tl]  using  (1.2)  and  the  current  estimates 
of  a  and  W  to  obtain  the  prior  bivariate  normal  distribution  of  (P,  \it) 

1  * 

having  mean  {bt_v  dMm M),  Var[p]  =  xt_v  Var[uf]  =  R(=    &t-\Ct_x  +  VVM, 

Corr(p,u,)  =    &t-iPt-iylci-i/Rr 

2.  Observe  the  Poisson  count  yt. 

3.  Invoke  the  Newton  procedure  of  Subsection  2.1  to  obtain  estimates  of 
the  parameters  of  the  posterior  bivariate  normal  distribution  of  (P,  \it) 

given  past  observations  and  the  new  observation  yf 

4.  Compute  new  estimates  of  a  and  W:   &t,Wr 

5.  Return  to  1. 

3.3  Results  of  Simulation  Experiments 

Various  simulation  experiments  were  carried  out  using  the  procedures 
outlined  in  this  section  to  estimate  a  and  W.  One  striking  phenomenon  was 
that  if  the  count  data  has  a  long  run  of  zeroes,  then  the  Kalman  procedure 
reacts  by  estimating  a  >  1  and  trying  to  estimate  \it  to  be  a  negative  number 
large  in  absolute  value.  This  tendency  caused  the  numerical  Newton 
procedure  discussed  in  Section  2  to  become  unstable,  and  rendered  the 
predicted  value  of  Xt  very  slow  to  respond  to  positive  counts  when  they  did 
eventually  occur.  This  behavior  was  ameliorated  by  arbitrarily  putting  a 
lower  bound  of  -2  on  flr  A  run  of  zero  counts  could  also  make  the  estimate  of 
W  become  close  to  zero.  A  small  value  of  W  makes  the  Kalman  filter 
unresponsive  to  count  changes.  This  behavior  was  mitigated  by  recomputing 


11 


estimates  of  a  and  W  every  5  =  10  time  units  rather  than  at  every  time,  and 
using  all  previous  times  t  as  in  (3.7)-(3.8).  The  estimates  of  a  and  W  were  also 
not  computed  if  the  last  10  observed  counts  were  zero.  Further  a  lower 
bound  of  0.1  was  arbitrarily  set  on  estimates  of  W .  Such  heuristics  are  not 
claimed  to  be  optimal,  but  are  needed  to  allow  the  procedures  to  behave 
suitably.  A  search  for  more  systematic  procedures  can  be  carried  out  in 
future  work. 

The  following  are  empirical  observations.  The  estimates  of  a  and  W  using 
p.{t)  =  m  t  tends  to  produce  practically  acceptable  estimates  of  a  but 
underestimates  W — it  tends  to  make  it  very  small.  This  behavior  is  not 
unexpected  since  m(can  be  thought  of  as  a  smoothed  estimate  of  the  \it  and  so 

will  have  a  smaller  variance  than  \ir    This  underestimation  of  VV  is  not 

11 

Neither  does  using  "windows"  of  length  8 


improved  by  using  fit  -  In 


Vt+2 


seem  to  solve  the  problem.    Other  procedures  for  estimating  a  and  W  are 
explored  in  the  following  two  sections. 

4.     ESTIMATION  OF  a  AND  W  FOR  THE  POISSON/NORMAL  MODEL: 
THE  FREEMAN-TUKEY  TRANSFORMATION 

This  section  reports  another  possible  approach  to  the  estimation  of  the 

autoregressive  parameters  for  a  Poisson/normal  model. 

The  model  is  as  before 

Hf+1  =  a\it  +cof+1  (4.1) 

where  {©,}  are  iid  normal  with  mean  0  and  variance  W; 

\e»t*frthAy 

P{%  =y^,xf}=exp{-^+Px'/2^ — i   .  (4.2) 

y ' 
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For  simplicity  of  notation  we  will  assume  ht  =  1. 

If  W  and  a  are  known,  then  an  approximate  Kalman  procedure  has  been 
given  previously.  The  issue  here  is  the  estimation  of  a  and  W. 

Let 

g{x)=\[x'  +  yfx~Tl  (4.3) 

and 

Z,  =  gOQ.  (4.4) 

The  conditional  distribution  of  Zt  given  \it  is  approximately  normal  with 

r     x      ,05 
mean  /:(p.(  +  x(P),  where  A:(x)  =  |4£    +  1J      ,  and  variance  1;  c.f.  Freeman  and 

Tukey  (1949, 1950)  and  Bishop,  Fienberg  and  Holland  (1975).  Thus  if  {Vf+1}  are 

iid  standard  normal,  given  Dt  =  (Y0,  Yv  ...,  Yf) 

d 
Z*+l"*(M*+l+*t+lP)+v*+l 

=  k{a\Lt  +  (&t+1  +  xt+iP)  +  Vf +1 

=  fc(aw%  +  a(\Lt  -mt)  +  (ot +1  +  x,+1fy  +  xf+i(P  -  bt))  +  V,+1 

*  fc(amt  +  *f+1fy)+fc'(amf  +  xt+1bt)[a(\it  -mt)+xt+1$-  bt)  +  (Ot+1\  +Vf+1.  (4.5) 

Hence, 

E[Z,+1|D,]  =  /;(am,  +  x,+1&,)  ,  (4.6) 

Var[Z,+1|D,]  =  [fc'(am(  +  xf+1fc,))2[a2Cf  +  xf+1Tt  +  2axt+1pt ^CtJ^t  +  W]  +1.  (4.7) 
An  approximate  joint  distribution  of  the  transformed  observations  is 


13 


P{Zi  £dzi,Z2  edz2,...,ZT+l  £dzT+i} 

=  P{Zl  ed21}P{Z2  £dz2\Zl  =  z1}*...  *P{ZT+1  Gd2T+1|Z!  =z1,...,ZT  =  zT} 

=  P{Y1  =  dy1}P{Z2  ed22|Y1=y1}x...xP{zT+1  e^^  =  yv...  ,YT  =  yT} 

-  p&i  =  yi}Il  ^=77  TAMo5exp|4[2'+i ' k{amt  +  ^+lfcf)'2  //*(a'w> I  <4-8> 
t«i  -J2nft(a, W)  I  ^  J 

where 

/f(a,W)  =[a2Q  +  x?+1Tf  +2ap^+1>/Q1/r7  +  w][/c'(amf  +  Xj+l^)]    +L 

Hence  an  approximate  In-likelihood  function  is  (up  to  addition  of  constants), 

T 
l{a,W)  =  £  --ln/£(a,W)  --(z^j  -  k{amt  +  xt+lbt))2ft(a,W)  m\  (4.9) 

t=\    2  2 


Note  that 


r      r        l05 

ifc(x)  =  [4ex  +  l]      ; 


0-5     r  2e 


X 


Jfc'(x)  =-[4ex+ll        4e*  =     ,  ; 

2L       J  viTTT 


2e*[2e*+l] 
[4e*  + 1] 


k"(x)  = ~—[^-  (4-!0) 


Further, 
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a  9 

ft(a,W)  =  [k'{amt  +  xt+1bt)\ 


dW 


(4.11) 


s  9 

— /,(a,W)  =  [2aQ  +  2xt+1ptJqJ^][k'(amt  +xt+1bt)] 

+[a2Ct  +x}+1Tt  +  2aPtxt+1JqJ^  +  W]2[k'{amt  +xt+1bt)] 


*k"(amt  +  xf+1frf)mf. 


(4.12) 


Differentiating  (4.9)  with  respect  to  a  and  W  we  obtain 

a 


dl      -     l^a 


ft(a,w) 


—  =  y  - 

aa     tfi    2    /t(«/W) 


1- 


i(«,VV) 


(4.13) 


1 

+(zt+l  -k(amt  +xt+1bt))k'{amt  +  xt+1bt)mtft{a,W)' 


dl 
dW 


1       1 


"£"  2/f(a,W)aW 


/f(a,W)  +  -[z,+i  -  fc(am(  +  *,+1&,)]  &*- 


ft(a,wy 


(4.14) 


=  V  -i —ft(a,W) 

£    2/f(a,W)d.W 


1  - 


[2f+1  -<:(afflt  +xt+ibt)\ 
ft{a,W) 


da2 


T 

El 


-a 


da 
/f(a,W) 


[k'(amt  +  xt+1bt)mt\ 
ft(a,W) 


(4.15) 


d2W 


T  1 


dW 


A(«/W) 


/f(a,W) 


(4.16) 
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■  a2/  " 

dadW  _ 

T  1 

aw 

f(a,W)^-/,(a,W) 

aa 

/,(a,W)2 

(4.17) 


To  avoid  difficulties  with  a  Newton  procedure  involving  the  restricted 
range  of  W  >  0,  we  will  reparameterize;  W  =  e.  In  this  case 

a/     a/   v 

—  = e7 


dy     dW 


Pr2] 


=  E 


1l 

dW7 


,2y 


\  d2l  1 
dady 

=  E 

a2/    1 

and  W  is  replaced  by  e  in  all  the  expressions  (4.12)  -  (4.17). 

A  Kalman  procedure  with  this  method  of  estimating  a  and  W  is  similar  to 
that  of  Section  3.2  except  that  Step  1  is  replaced  by  the  numerical  solution  of 
the  system  of  equations, 


da  =  ° 


d-y  =  ° 


(4.18) 


(4.19) 


by  a  Newton  procedure  using  Fisher's  scoring  method.  Set  W  =  e r  where  y  is 

the  solution.   Occasionally,  there  are  numerical  problems  with  the  Newton 

procedure.   In  these  cases  search  is  used  to  find  estimates  of  a  and  W. 

A  summary  of  the  Newton  (with  default  search)  procedure  is  as  follows. 

The  procedure  starts  with  the  previously  estimated  a  and  /called  a0  and  yQ. 
0.     Set  B  =  0. 
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1.     Compute    \f(S/,ea)  ,  \f(a/,ay)  ,  £  \b\bc\[(\£(e2!,da2)),  £ 


r4 

dy 
dy 


■    and  E[dady    using  ao  and  ?V    If 
<  O.OOOl  go  to  2',  otherwise  go  to  2. 


'-£-1 

_dady_ 


<  0.001  or  E 


2.     Solve  the  system  of  linear  equations  (in  terms  a  and  y) 


6a 


da2 


(a  -a0)  +  E 


eV 

dady 


(r  -y0) 


67 


dady 


a  -  a0)  +  E 


dy2 


(r  -  r0) 


for  an  and  yn.  Go  to  3. 

r)/ 

2'.    Set  yn  =  yQ.  If  x~  is  of  the  same  sign  for  both  the  current  and  previous 
value  of  a,  set 


a/ 

an-aQ-—l 


da2 


■at 

and  go  to  3.    If  ^t  is  of  different  signs  for  the  current  and  previous 
values  for  a,  then  a  golden  section  search  is  used  to  solve  the  equation 

da 


set  an  equal  to  the  first  value  of  a  for  which 


dl 


da 


<  0.1;  if  the  number  of 


iterations  for  this  one-dimensional  search  is  greater  than  50,  go  to  4. 


3.     If  max 
values 


dl 


da 


(<*J, 


dl 


da 


(r») 


<  0.1  stop  and  take  an  and  yn  as  the  estimated 


3'.    If  I  an  I  <  5  and  yn  <  5,  go  to  5. 
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3".    If  |  a  J  >  5  or  yn  >  5,  let  B  =  B+l.  If  B  =  1,  go  to  4.   If  B  =  2,  go  to  4'.  If 

B  >  3,  to  to  4". 

4.  Do  a  two-dimensional  search  of  the  ln-likelihood  function  (4.9) 
parameterized  in  terms  of  a  and  W  for  a  maximizing  value.  The  grid 
for  a  is  [-4,4]  in  steps  of  0.2;  the  grid  for  W  is  [0,8]  in  steps  of  0.2.  If 


max 


dl 


da 


(««)/ 


dl 


dy 


(InW) 


<0.1 


return  an  and  yn  =  InW  as  the  estimated  values;  otherwise  set  a0  =  an 
and  y0  =  InW  and  return  to  1. 

4'.  Same  as  4  except  the  grid  for  a  is  [-8,8]  in  steps  of  0.1  and  the  grid  for  W 
is  [0.1,8]  in  steps  of  0.1  with  the  additional  points  0.000001,  0.00001, 
0.0001,  0.001,  0.01. 

4".    y  is  set  equal  to  its  previous  value  and  a  golden  section  search  for  the 

r  dl    ■ 

zero  of  t~  is  done  as  in  2  . 

5.     Set  ao  -  ctn  and  y0  =  yn  and  return  to  1  if  the  number  of  iterations  is  less 

than  50.  If  the  number  of  iterations  is  greater  than  50  set  B  =  B  +  1.   If 
B  =  1,  go  to  4.    If  B  >  6,  return  an  and  yn  as  the  estimated  values.    If 

1  <  B  <  5,  go  to  6. 


dl 


and 

-<y> 
ey17"' 

• If  1 

do^K) 

<     | 

~(y) 

go  to  6'; 


6.     Compute    ^(«„) 
otherwise  go  to  6". 
6'.    Fix  an  -  aQ  and  do  a  search  over  the  interval  [-4B,  4B]  for  that  yfor 


which 


$1 
dy' 


dl 


<  0.1;  if  the  number  of 


Set  yn  equal  to  the  first  value  of  y  for  which  L^ 
iterations  is  greater  than  50,  go  to  4. 
6".   Fix  yn  =  y0  and  do  a  search  over  the  interval  [-4B,  4B]  for  that  a  for 


dl 


which  0  =  -jr~ ;  set  an  equal  to  the  first  value  of  a  for  which 
if  the  number  is  greater  than  50,  go  to  4. 


dl 


da 


<0.1; 
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The  following  are  empirical  observations.  The  procedure  involves  a  great 
deal  of  computing.  For  small  times  t,  the  estimate  of  W  can  be  very  close  to 
zero;  a  lower  bound  of  0.1  was  placed  on  the  value  of  W  that  is  input  into  the 
Kalman  procedure.   Further,  a  lower  bound  of  -2  was  also  put  on  p...   The 


stopping  criteria  of  max 


dl 


da 


— 
'  dy 


<  0.1  seems  to  be  adequate;  a  smaller 
tolerance  can  result  in  the  procedure  becoming  unstable  and  greatly 
increasing  the  computational  effort. 

5.     ESTIMATION  OF  a  AND  W  FOR  THE  POISSON/NORMAL  MODEL: 
LOGARITHMIC  TRANSFORMATION 

Results  of  Lambert  (1989)  suggest  that  when  count  data  arise  from  a 
mixture  of  Poisson  distributions,  then  a  smaller  power  transformation  of  the 
counts  than  the  square  root  is  needed  to  stabilize  the  variance  of  the  count 
data.  In  particular,  if  the  variability  of  the  mixture  is  largish,  then  a  log 
transformation  may  stabilize  the  variance.  In  this  section  simple  procedures 
for  estimating  a  and  W  based  on  a  log  transformation  of  the  count  data  are 
described. 

Suppose  \t  and  \it  are  as  in  (1.1)  and  (1.2).  Given  \it+v  (J,  the  first  two  terms 
of  a  Taylor  expansion  yield 

lnY, +1  *  lnfc, +1  +  % +1p  +  jif  +1  +  [htexp{xt+1$  +  nt+1}] '  [Y,  -  ktexp{xt+1$  +  ^+1}] 

(5.1) 

Let 


Z,+1  =ln 


Vl+" 


-  lnht+1  -  xt+i$. 


(5.2) 


Using  the  first  term  of  the  Taylor  expansion  (5.1) 
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zf+l  *  m+l  =  a\k  +  ©t+l  "  amt  +  mt+l-  (5.3) 

Hence,  given  [Yq,  Yj,  ...,  Yf],  the  distribution  of  Z/+1  has  approximate  mean  am, 

and  variance  W.    We  will  approximate  the  distribution  of  Zt+1  with  a  normal 

distribution  having  mean  am t  and  variance  W.  A  ln-likelihood  function 

under  this  approximation  is  (up  to  addition  of  constants) 

T  -1 
l(a,W)  =  £  -±lnJV  .i(2fc+1  -  am,)2-!-;  (5.4) 

t-l     2  2  w 

a/       T - 1  1 

7-  =  E  (2'+i  "  am*  K  T77  =  0;  (5>5) 

ai'  I;1  n    i  ,2  i 

=  V +  -(zt+1  -  amt)  — T  =  0.  (5.6 

ew      £i    2VV      2V  '  *  f     VV2 

Thus, 

r-i 

X  (Zf+1%) 

«Lm=-^TTi (5-7) 

5>2 
t-i 

1     T-!  2 

WL(T)=  — -  £  (zf+1  -  dL(T)mf)  (5.8) 

are  the  resulting  estimates  of  a  and  W. 

Simulation  experiments  suggest  that  dL  and  W^  tend  to  have  a  large 
positive  bias. 

Another  possibility  is  to  estimate  a  by 

-T-l  T-l 

«m(T)  =  7  E  ™t  +1™,  /  I  m?  (5'9) 

1    t-l  t-l 

as  in  (3.7)  of  Section  3  and  W  by 
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w 


1  T-1 


^T)  =  777  I  (2*+l  -  OmCD^)  •  (5-1Q) 


m 

These  estimates  are  easy  to  compute.  Simulation  experiments  suggest  that 
this  estimate  of  am  tends  to  have  a  negative  bias  and  Wm  L  has  a  positive  bias 
but  not  as  large  as  that  of  VV^ 

6.     RESULTS  FROM  SIMULATION  EXPERIMENTS 

In  this  section  we  report  results  of  simulation  experiments  concerning  the 
behavior  of  the  approximate  Kalman  filter  of  Section  2  under  various 
procedures  for  estimating  a  and  W. 

Each  replication  of  the  simulation  consists  of  generating  a  Poisson  time 
series  {yt;  t  -0,1,  2, ...,  T}  using  equations  (1.1)  -  (1.2).  The  random  numbers 
were  generated  using  LLRANDOM  II,  cf.  Lewis  et  al.  (1981).  If  a  <  1,  then  jig 
is  drawn  from  the  stationary  distribution  of  {\it;  t  >  0}.  The  Kalman  procedure 
with  each  of  the  procedures  for  estimating  a  and  W  is  then  applied  to 
{y{  t  =  0,1,  2, ...,  T}.  For  each  time  t,  the  estimate 

Xt  =  exp\htexp{{Lt  +  xj}}  (6.1) 

is  computed  and  the  mean  square  prediction  error 

T"1  .2 


1  ' L  t=i 


is  calculated.  In  addition  the  average  estimated  values  of  a  and  W  and  /3  are 
computed: 

1  T 
m(a)  =  -£a,  (6.3) 

1  *=l 
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1  T 

m(W)=-£w,  (6.4) 

1  t=l 


1  T 
m(/3)  =  -J>  (6.5) 


where  dt  and  VVf  are  the  estimates  of  a  and  W  obtained  after  the  observation  at 
time  t  and  bt  is  as  in  Section  1. 

The  simulation  is  replicated  N  times.  If  e^m^a),  m  ,-(W),  m  -(/J)  represent 
the  summary  statistics  from  the  *  replication,  then  the  mean  of  the  summary 
statistics  are  computed  for  the  N  replications;  i.e. 

1  N 
e=  —  j^ep  (6.5) 


N 

I 

N 


I 1       N  2 


a  =-i-j;m/(a);  (6.7) 

N  i=\ 

-       1   N 

w"77l>i(W);  (6-8) 

N  i=\ 

-       1    N 

/?=-5>,(/3).  (6.8) 

These  latter  statistics  are  then  compared  for  the  different  procedures  of 

estimating  a  and  W. 

Results  for  the  following  procedures  for  estimating  a  and  W  are  reported. 
1.     DD:   a  and  W  are  both  known  and  are  not  estimated. 
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2.  FT:  Both  a  and  W  are  estimated  by  a  two-dimensional  search  of  the  ln- 
likelihood  based  on  the  Freeman-Tukey  transformation  (4.9).  The  grid 
for  W  is  [0.1,3]  in  increments  of  0.1;  the  grid  for  a  is  [-3,3]  in  increments 
of  0.1. 

3.  M/FT:  a  is  estimated  by  the  moment  estimator  using  (3.7);W  is 
estimated  by  searching  the  one-dimensional  likelihood  based  on  the 
Freeman-Tukey  transformation  with  &  set  equal  to  (3.7);  the  grid  is 
[0.1,3]  with  steps  of  size  0.1. 

4.  M/LN:  a  is  estimated  using  the  moment  estimate  of  (3.7).  Wis 
estimated  using  the  moment  estimator  (5.10)  based  on  the  logarithmic 
transformation. 

5.  LN/LN:  a  and  W  are  estimated  using  the  moment  estimators  (5.9)  and 
(5.10)  based  on  the  logarithmic  transformation. 

All  the  procedures  have  a  lower  bound  of  -2  on  fit,  a  lower  bound  of  0.1 
on  VVp  and  an  upper  bound  of  1  on  the  absolute  value  of  &t. 

For  the  results  of  the  simulation  experiments  reported  in  Table  1,  a  =  0.5, 
W  =  0.25,  xt*l,ht*l,fi  =  0.5.  The  initial  values  of  the  estimates  VV0  =  0.25,  &6  = 

0.5,  fi0  =  0;  for  the  Kalman  C0  =  1,  t0  =  1,  /50  =  0,  $  =  0.   The  two-dimensional 

search  for  the  estimates  of  a  and  W  in  FT  takes  a  large  amount  of  computer 

time.  As  a  result,  the  number  of  replications  for  the  experiments  reported  in 

Table  1  is  small.   However,  the  results  of  Table  1  suggest  that  the  two  most 

promising  procedures  are  M/LN  and  M/FT. 

In  Table  2,  the  number  of  replications  is  100.  The  procedures  M/LN  and 
M/FT  only  are  compared.  In  Table  2,  the  initial  values  of  VV0  =  0.25,  &0  =  0.5, 

fi0  =  0;  for  the  Kalman  C0  =  1,  %  =  1,  fa  =  0,  fi  =  0.  . 

In  Table  3,  {xt}  =  (0.25,  0.5, 1,  0.25,  0.5, 1,  ...},  j3  =  0.5,  W  =  0.25,  a  =  0.5.  The 
initial  values  of  VV0  =  0.25,  &0  =  0.5,  fiQ  =  0;  for  the  Kalman  C0  =  1,  t0  =  1,  /50  =  0, 

0  =  0. 

In  Table  4,  the  parameters  are  the  same  as  in  Table  3  except  W  =  1.  The 
initial  values  of  VV0  =  0.25,  &0  =  0.5,  fi0  -  0;  for  the  Kalman  C0  =  1,  f0  =  1,  /50  =  0, 

0  =  0. 
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Comparison  of  Tables  2  and  3  suggests,  not  surprisingly,  that  the  mean 
square  prediction  error  is  smaller  when  there  is  a  variable  covariate  {*,}. 

Comparison  of  Tables  3  and  4  suggests,  not  surprisingly,  that  a  larger 
value  of  W  results  in  a  larger  mean  square  prediction  error. 
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TABLE  1. 


SERIES  LENGTH  =  5;  20  REPLICATIONS 


METHOD: 

DD 

FT 

M/FT 

M/LN 

LN/LN 

Mean  MSE 

4.34 

4.91 

4.77 

4.19 

4.96 

St.  Dev.  MSE 

3.21 

3.82 

3.64 

3.17 

3.76 

Max  MSE 

11.0 

13.2 

11.7 

12.5 

13.5 

Mean  6c 

- 

0.10 

0.11 

0.15 

0.26 

Mean  W 



0.48 

0.67 

0.98 

0.77 

SERIES  LENGTH  =  10;  20  REPLICATIONS 


METHOD: 

DD 

FT 

M/FT 

M/LN 

LN/LN 

Mean  MSE 

2.94 

3.15 

3.08 

2.77 

3.36 

St.  Dev.  MSE 

2.27 

2.42 

2.28 

1.93 

2.54 

Max  MSE 

10:4 

11.1 

10.1 

8.27 

11.2 

Mean  6c 

- 

0.12 

0.08 

0.09 

0.19 

Mean  VV 



0.23 

0.32 

0.78 

0.69 

SERIES  LENGTH  =  20;  40  REPLICATIONS 


METHOD: 

DD 

FT 

M/FT 

M/LN 

LN/LN 

Mean  MSE 

3.38 

3.90 

3.72 

3.29 

3.97 

St.  Dev.  MSE 

2.00 

2.81 

2.43 

2.16 

2.62 

Max  MSE 

10.7 

13.5 

12.1 

10.3 

11.2 

Mean  6c 

- 

0.33 

0.15 

0.17 

0.31 

Mean  VV 



0.22 

0.31 

0.79 

0.73 
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TABLE  2.  xtszl,a  =  0.5,  W  =  0.25 
SERIES  LENGTH  =  5;  100  REPLICATIONS 


METHOD: 

DD 

M/FT 

M/LN 

Mean  MSE 

4.21 

4.75 

4.32 

St.  Dev.  MSE 

4.68 

5.81 

5.40 

Max  MSE 

28.7 

38.0 

34.1 

Mean  & 

- 

0.08 

0.09 

Mean  VV 

- 

0.55 

0.93 

Mean  /3 

0.34 

0.25 

0.23 

SERIES  LENGTH  =  10;  100  REPLICATIONS 


METHOD: 

DD 

M/FT 

M/LN 

Mean  MSE 

3.98 

4.41 

4.14 

St.  Dev.  MSE 

3.27 

3.87 

3.58 

Max  MSE 

16.4 

19.0 

17.3 

Mean  & 

- 

0.07 

0.10 

Mean  VV 

.• 

0.45 

0.88 

Mean  /3 

0.42 

0.36 

0.31 

SERIES  LENGTH  =  20;   100  REPLICATIONS 


METHOD: 

DD 

M/FT 

M/LN 

Mean  MSE 

3.96 

4.35 

4.18 

St.  Dev.  MSE 

2.59 

3.30 

3.06 

Max  MSE 

20.3 

26.5 

23.6 

Mean  & 

- 

0.10 

0.13 

Mean  VV 

- 

0.41 

0.82 

Mean  (i 

0.48 

0.43 

0.37 
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TABLE  3.  VARIABLE  {*,},  a  =  0.5,  W  =  0.25 


SERIES  LENGTH  =  5; 

100  REPLICATIONS 

METHOD: 

DD 

M/FT 

M/LN 

Mean  MSE 

3.69 

3.89 

3.70 

St.  Dev.  MSE 

3.83 

4.09 

4.00 

Max  MSE 

15.6 

18.4 

19.9 

Mean  & 

- 

0.13 

0.12 

Mean  VV 

- 

0.52 

0.80 

Mean  J3 

0.32 

0.21 

0.18 

METHOD: 

Mean  MSE 
St.  Dev.  MSE 
Max  MSE 
Mean  & 
Mean  VV 
Mean  /3 


SERIES  LENGTH  =  20;  100  REPLICATIONS 

DD  M/FT  M/LN 

3.11  3.32                             3.21 

2.48  2.80                             2.68 

20.4  22.7  20.2 

0.12  0.10 

,.       -  0.35                             0.78 

0.41  0.36                             0.28 
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TABLE  4.   VARIABLE  {*,},  a  =  0.5,  W  =  1 
SERIES  LENGTH  =  5;   100  REPLICATIONS 


METHOD: 

DD 

M/FT 

M/LN 

Mean  MSE 

19.21 

20.70 

19.88 

St.  Dev.  MSE 

52.2 

53.0 

52.1 

Max  MSE 

415.1 

390.1 

393.5 

Mean  & 

- 

0.19 

0.18 

Mean  VV 

- 

0.93 

1.05 

Mean  (5 

0.31 

0.14 

0.12 

SERIES  LENGTH  =  20;  100  REPLICATIONS 


METHOD: 

DD 

M/FT 

M/LN 

Mean  MSE 

18.93 

21.6 

20.4 

St.  Dev.  MSE 

28.6 

31.5 

30.9 

Max  MSE 

181.7 

192.1 

188.5 

Mean  & 

- 

0.20 

0.19 

Mean  VV 

- 

1.2 

1.1 

Mean  /3 

0.55 

0.32 

0.31 
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Of  the  two  procedures,  M/LN  tends  to  have  the  smaller  mean  MSE. 
M/LN  tends  to  have  a  positive  bias  estimating  W  and  negative  biases 
estimating  &  and  /3.  The  procedure  M/FT  takes  more  computational  effort, 
tends  to  underestimate  W,  a  and  /?  and  produces  slightly  higher  mean  MSE. 

Tables  5-7  report  results  comparing  procedures  estimating  a  and  W  with 

the  procedure  NFT  in  which  both  a  and  W  are  estimated  using  the  Newton 

procedure  of  Section  4.  For  all  the  procedures  reported  in  these  tables,  there 

is  no  bounding  on  &.   In  Table  5,  xt  ■  1,  ht  ■  1,  f$  =  0.5,  a  =  0.5,  W  =  0.25.  In 

Table  6,  xt  =  0.25,  0.5, 1,  0.25,  0.5, 1, ...,  and  the  other  parameters  are  as  before. 

For  Tables  5  and  6  the  initial  values  of  VV0  =  0.25,  d0  =  0.5,  fi0  =  0,  C0  =  1,  f0  =  1, 

pQ  =  0.  The  values  in  parentheses  are  estimates  of  standard  deviations;  e.g. 

1 
1      ",      .   .     -,2~ 


—Yintila)  -a)' 


For  Table  7,  {rf}  is  as  In  Table  6,  a  =  0.5,  W  =  0.25,  f$0  =  0.5,  and  the  initial 
estimates  are  &0  =  0,  VV0  =  1,  C0  =  1,  t0  =  1,  pQ  =  0,  /20  =  0. 

The  Newton  procedure  of  Section  4  does  not  take  as  much  time  as  the 
two-dimensional  search  used  in  Tables  1-2;  however,  it  is  still  a  much  larger 
computational  effort  than  the  other  two  procedures.  Once  again,  based  on  the 
mean  MSE  and  computational  effort  the  procedure  M/LN  appears  to  be  the 
most  attractive.  The  procedure  M/LN  once  again  appears  to  overestimate  W 
and  underestimate  a  and  /3. 
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TABLE  5.  xt  s  1;  a  =  0.5,  W  =  0.25,  0  =  0.5 


SERIES  LENGTH  =  5,  25  REPLICATIONS 


METHOD: 

DD 

NFT 

M/FT 

M/LN 

Mean  MSE 

3.94 

4.44 

4.32 

3.89 

Std  Dev.  MSE 

3.08 

3.66 

3.49 

3.10 

Max  MSE    . 

11.0 

13.2 

11.7 

12.5 

Mean  & 

- 

0.06  (0.37) 

0.08  (0.22) 

0.10  (0.23) 

Mean  VV 

- 

0.46  (0.47) 

0.60  (0.60) 

0.93  (0.45) 

Mean  $ 

0.28 

0.18 

0.21 

0.18 

SERIES  LENGTH  =  20,  25  REPLICATIONS 

METHOD: 

DD 

NFT 

M/FT 

M/LN 

Mean  MSE 

3.11 

3.40 

3.32 

3.14 

Std  Dev.  MSE 

1.83 

2.12 

2.11 

2.03 

Max  MSE 

8.50 

9.29 

9.60 

8.78 

Mean  a 

- 

0.22  (0.39) 

0.10  (0.30) 

0.12  (0.29) 

Mean  W 

- 

0.29  (0.24) 

0.32  (0.26) 

0.80  (0.23) 

Mean  J3 

0.38 

0.23 

0.31 

0.23 
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TABLE  6.  VARIABLE  {*,} ;  a  =  0.5,  W  =  0.25,  0  =  0.5 


SERIES  LENGTH  =  20,  25  REPLICATIONS 

METHOD: 

DD 

NFT 

M/FT 

M/LN 

Mean  MSE 

2.63 

2.77 

2.81 

2.63 

Std  Dev.  MSE 

1.78 

1.95 

2.06 

1.79 

Max  MSE 

9.24 

10.10 

10.86 

9.12 

Mean  & 

- 

0.21  (0.38) 

0.06  (0.32) 

0.06  (0.34) 

Mean  VV 

- 

0.20  (0.14) 

0.25  (0.20) 

0.68  (0.18) 

Mean  /} 

0.33 

0.21 

0.27 

0.18 
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TABLE  7.   VARIABLE  {*,} ;  a  =  0.5,  W  =  0.25,  0  =  0.5 
d0  =  0,W0  =  l/C0=l,  i0  =  l,,50  =  0,  00  =  0 

SERIES  LENGTH  =  20, 100  REPLICATIONS 


METHOD: 

DD 

NFT 

M/FT 

M/LN 

Mean  MSE 

2.95 

3.25 

3.21 

3.04 

Std  Dev.  MSE 

2.21 

2.70 

2.67 

2.28 

Max  MSE 

13.2 

16.5 

16.5 

13.8 

Mean  & 

- 

0.12  (0.34) 

0.11  (0.28) 

0.09  (0.29) 

Mean  W 

- 

0.34  (0.44) 

0.36  (0.32) 

0.73  (0.19) 

Mean  $ 

0.36 

0.25 

0.28 

0.21 
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7.     SUMMARY  AND  CO  VCLUSIONS 

In  this  paper  we  have  investigated  approximate  methods  for  estimation 
and  prediction  for  the  hierarchical  Poisson/normal  time  series  model  given 
in(l.l)-(1.2).  For  given  values  of  the  random  walk  parameters,  a  and  W,  the 
joint  distribution  of  (P,  uf)  is  approximated  by  a  bivariate  normal  distribution 
using  the  Laplace  method.  Various  heuristic  methods  for  estimating  a  and  W 
are  presented.  Based  on  simulation  results  and  ease  of  computation,  it  is 
recommended  that  a  be  estimated  by  (3.7) 

T  T 

a(T)  =  Y,mtmt-l/Y,m?-i 

and  W  be  estimated  by  (5.10) 

1     TJ,  ^       ,2 


wm=  — £(z,+1-a(T)m,r 


where  zt+1  is  given  by  (5.2). 

These  estimates  of  a  and  W  along  with  the  approximate  Kalman 
procedure  using  the  Laplace  method  provide  a  computationally  easy 
procedure  for  prediction  of  the  time  series. 

Figures  4-7  show  results  of  a  simulation  of  model  (1.1)-(1.2)  for 
1  =  1,.,.,  100.  In  the  example  /3  =  0.5,  {xt}  takes  the  values  {0.25,  0.5, 1}  over  and 
over,  and  ht  =  1.  {(ot}  are  iid  normal  with  mean  0,  variance  0.25,  and  a  =  0.5. 

The  simulation  starts  with  \Iq  drawn  from  the  stationary  distribution  of  {\Lt},  a 

2       1 
normal  distribution  with  mean  0  and  variance  W/(l-a  )  =  ~.  The  values  of  a 

and  W  are  estimated  using  (3.7)  and  (5.10).  Initial  values  of  a  and  W  are 
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a0  =  0.5  and  VV  =  1,  C0  =  1,  t0  =  1,  p0  =  0,  /30=  0,  /x0  =  0.  There  are  no  bounds  on 
flt,  a  and  VV  . 

Figure  4  displays  the  count  data.  Also  displayed  are  the  true  Xt  (circles), 
the  estimated  Xt  when  both  a  and  VV  are  known  (dashed  line)  and  the 
estimated  Xt  when  both  a  and  VV  are  estimated  (solid  line).  The  Xt  when  both  a 
and  VV  are  estimated  is  perhaps  a  little  more  responsive  to  changes  in  the  data 
than  when  a  and  VV  are  known.  The  difference  between  the  two  estimates  of  Xt 
is  greatest  for  small  times  t. 

Figure  5  presents  plots  of  the  estimated  Vw  and  a.  Note  that  VV  has  a 
positive  bias  which  will  make  the  Kalman  procedure  more  responsive  to  the 
count  data.  The  estimated  values  of  a  appear  to  be  reasonable. 

Figure  6  presents  estimates  of  j3  and  its  estimated  standard  deviation  yfc 
both  for  the  Kalman  with  a  and  VV  known  (dotted  line)  and  with  a  and  VV 
unknown  (solid  line).  The  estimates  of  Tf  generally  decrease  as  t  increases 
reflecting  the  model  assumption  that  /3  is  constant.  The  estimates  of  fi  appear 
to  be  reasonable.  The  estimates  of  /3  with  a  and  VV  also  estimated  tend  to  be 
smaller  than  those  for  which  a  and  VV  are  known;  this  behavior  may  be  due  to 
the  fact  that  estimation  of  a  and  VV  is  accounting  for  some  of  the  variability  of 
the  data  that  would  otherwise  be  accounted  for  by  /3. 

Figure  7  displays  the  true  \it  (circles),  the  estimated  ]i t  =  m  t  with 
parameters  a  and  W  known  (dotted  line)  and  the  estimated  fit  with  parameters 
a  and  W  estimated  (solid  line).  The  count  data  are  also  displayed.  The 
estimates  of  fit  with  a  and  W  estimated  are  more  variable  than  those  when  a 
and  W  are  known  reflecting  the  greater  responsiveness  of  the  Kalman  to  the 
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data  when  the  estimated  W  has  a  position  bias.  Once  again  the  estimates  of  pt 
appear  to  be  reasonable. 
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